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Abstract: In this article we propose a shooting algorithm for a class of opti- 
mal control problems for which all control variables appear linearly. The shoot- 
ing system has, in the general case, more equations than unknowns and the 
Gauss-Newton method is used to compute a zero of the shooting function. This 
shooting algorithm is locally quadratically convergent if the derivative of the 
shooting function is one-to-one at the solution. The main result of this paper is 
to show that the latter holds whenever a sufficient condition for weak optimality 
is satisfied. We note that this condition is very close to a second order necessary 
condition. For the case when the shooting system can be reduced to one having 
the same number of unknowns and equations (square system) we prove that the 
mentioned sufficient condition guarantees the stability of the optimal solution 
under small perturbations and the invertibility of the Jacobian matrix of the 
shooting function associated to the perturbed problem. We present numerical 
tests that validate our method. 

Key-words: optimal control, Pontryagin Maximum Principle, singular con- 
trol, constrained control, shooting algorithm, second order optimality condition, 
stability 
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Un algorithme de tir pour les problemes de 
commande optimale avec des arcs singuliers 

Resume : Dans ce travail on presente une condition sufhsante pour que 
l'algorithme de tir soit localement convergent quand il est applique aux problemes 
de commande optimale afnnes dans les commandes. On commence par etudier le 
cas avec des contraintes initiales-finales sur l'etat et commande librc, ct en suite 
on ajoutc des contraintes sur la commande. L'algorithme de tir est localement 
quadratiquement convergent si la derivee de la fonction de tir associee est 
injective dans la solution optimale. Le resultat principal de cet article montre 
une condition sufnsante pour cette injectivite, qui est tres proche de la condition 
necessairc du second ordre. On montre que cette condition sufnsante assure la 
stabilite de la solution optimale aux petites perturbations ct qu'elle garantit 
aussi que l'algorithme de tir est convergent pour le probleme perturbe. On 
presente des essais numcriqucs qui validcnt notre methode. 

Mots-cles : commande optimale, Principe de Pontryaguinc, commande 
singulicrc, contraintes sur la commande, algorithme de tir, conditions d'optimalitc 
du second ordre, stabilite 
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1 Introduction 

The classical shooting method is used to solve boundary value problems. Hence, 
it is used to compute the solution of optimal control problems by solving the 
boundary value problem derived from the Pontryagin Maximum Principle. 

Some references can be mentioned regarding the shooting method. The 
first two works we can find in the literature, dating from years 1956 and 1962 
respectively, are Goodman-Lance [1] and Morrison et al. [2]. Both present the 
same method for solving two-point boundary value problems in a general setting, 
not necessarily related to an optimal control problem. The latter article applies 
to more general formulations. The method was studied in detail in Keller's 
book [3J , and later on Bulirsch [J] applied it to the resolution of optimal control 
problems. 

The case we deal with in this paper where the shooting method is used to 
solve optimal control problems with control-affine systems is treated in, e.g., 
Maurer [5], Oberle [DJ [7], Fraser- Andrews [5], Martinon [pj and Vossen [ID] . 
These works provided a series of algorithms and numerical examples with dif- 
ferent control structures, but no theoretical foundation is supplied. In particu- 
lar, Vossen [10] dealt with a problem in which the control can be written as a 
function of the state variable, i.e. the control has a feedback representation. He 
proposed an algorithm that involved a finite dimensional optimization problem 
induced by the switching times. The main difference between Vossen's work and 
the study here presented is that we treat the general problem (no feedback law 
is necessary). Furthermore we justify the well-definition and the convergence 
of our algorithm via second order sufficient conditions of the original control 
problem. In some of the just mentioned papers the control variable had only 
some of its components entering linearly. This particular structure is studied 
in more detailed in Aronna [11] , and in the present article we study problems 
having all affine inputs. 

In [12] Bonnard and Kupka studied the optimal time problem of a generic 
single-input affine system without control constraints, with fixed initial point 
and terminal point constrained to a given manifold. For this class of problems 
they established a link between the injectivity of the shooting function and the 
optimality of the trajectory by means of the conjugate and focal points theory. 
Bonnard et al. [13] provides a survey on a series of algorithms for the numerical 
computation of these points, which can be employed to test the injectivity of the 
shooting function in some cases. The reader is referred to [13], Bonnard-Chyba 
[14] and references therein for further information about this topic. 

In addition, Malanowski-Maurer [IS] and Bonnans-Hermant [ID] dealt with 
a problem having mixed control-state and pure state running constraints and 
satisfying the strong Legendre-Clebsch condition (which is not verified in our 
affine- input case). They all established a link between the invertibility of the 
Jacobian of the shooting function and some second order sufficient condition for 
optimality. They provided stability analysis as well. 

We start this article by presenting an optimal control problem affine in the 
control, with terminal constraints and free control variables. For this kind of 
problem we state a set of optimality conditions which is equivalent to the Pon- 
tryagin Maximum Principle. Afterwards, the second order strengthened gener- 
alized Legendre-Clebsch condition is used to eliminate the control variable from 
the stationarity condition. The resulting set of conditions turns out to be a two- 
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point boundary value problem, i.e. a system of ordinary differential equations 
having boundary conditions both in the initial and final times. We define the 
shooting function as the mapping that assigns to each estimate of the initial val- 
ues, the value of the final condition of the corresponding solution. The shooting 
algorithm consists of approximating a zero of this function. In other words, the 
method finds suitable initial values for which the corresponding solution of the 
differential equation system satisfies the final conditions. 

Since the number of equations happens to be, in general, greater than the 
number of unknowns, the Gauss-Newton method is a suitable approach for solv- 
ing this overdetermined system of equations. The reader is referred to Dennis 
[T7j . Fletcher [TS] and Dennis et al. [TH] for details and implementations of 
Gauss-Newton technique. This method is applicable when the derivative of the 
shooting function is one-to-one at the solution, and in this case it converges 
locally quadratically. 

The main result of this paper is to provide a sufficient condition for the in- 
jectivity of this derivative, and to notice that this condition is quite weak since, 
for qualified problems, it characterizes quadratic growth in the weak sense (see 
Dmitruk [2Ql EI])- Once the unconstrained case is investigated, we pass to a 
problem having bounded controls. To treat this case, we perform a transforma- 
tion yielding a new problem without bounds, we prove that an optimal solution 
of the original problem is also optimal for the transformed one and we apply 
our above-mentioned result to this modified formulation. 

It is interesting to mention that by means of the latter result we can justify, in 
particular, the invertibility of the Jacobian of the shooting function proposed by 
Maurer [5] . In this work, Maurer suggested a method to treat problems having 
scalar bang-singular-bang solutions and provided a square system of equations 
(i.e. a system having as many equations as unknowns) meant to be solved by 
Newton's algorithm. However, the systems that can be encountered in practice 
may not be square and hence our approach is suitable. 

We provide a deeper analysis in the case when the shooting system can 
be reduced to one having equal number of equations and unknowns. In this 
framework, we investigate the stability of the optimal solution. It is shown 
that the above-mentioned sufficient condition guarantees the stability of the 
optimal solution under small perturbation of the data and the invertibility of 
the Jacobian of the shooting function associated to the perturbed problem. 
Felgenhauer in H3] provided sufficient conditions for the stability of the 
structure of the optimal control, but assuming that the perturbed problem had 
an optimal solution. 

Our article is organized as follows. In section [2] we present the optimal 
control problem without bound constraints, for which we provide an optimality 
system in section [3j We give a description of the shooting method in section |4j 
In section[5]we present a set of second order necessary and sufficient conditions, 
and the statement of the main result. We introduce a linear quadratic optimal 
control problem in section [6j In section [7] we present a variable transformation 
relating the shooting system and the optimality system of the linear quadratic 
problem mentioned above. In section [8] we deal with the control constrained 
case. A stability analysis for both unconstrained and constrained control cases 



is provided in section Finally we present some numerical tests in section 10 
and we devote sectionTTl] to the conclusions of the article. 
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2 Statement of the Problem 

Consider the spaces U := £«,((), T; R m ) and X := W^(0, T; R"), as control and 
state spaces, respectively. Denote by u and x their elements, respectively. When 
needed, put w = (x, u) for a point in the product space W :— X x U. In this 
paper we investigate the optimal control problem 

J := ip (x , x T ) -> min, (1) 

m 

it = ^Ui, t /i(£ t ), a.e. on [0,T], (2) 

i=0 

rjj(x ,x T ) = 0, for j = 1, ...,<i n , (3) 

where final time T is fixed, uq = 1, /j : R n — > R™ for i — 0, . . . ,m and 7^ : 
R 2 ™ — > R for j = 1, ...,d v . Assume that data functions ipo, fi and rjj have 
Lipschitz-continuous second derivatives. Denote by (P) the problem defined by 
Q-Q. An element w € W satisfying ([2])-(|3]) is called a feasible trajectory. 

Set X* := 14^(0, T; R n '*) the space of Lipschitz-continuous functions with 
values in the n— dimensional space of row- vectors with real components M™'*. 
Consider an element A := (f3,p) € R dr "* x X* and define the pre-Hamiltonian 
function 

m 

F[A](a;,u,t) :=pt^Wi/i(i), (4) 

i=0 

the initial-final Lagrangian function 

*[A](Co,Ct) := ^o(Co, Ct) + £>»&(G>,Ct), (5) 

i=l 

and the Lagrangian function 

C[X](w):=£[X\(xo,x T )+ I Pt I y ]ui,tfi(x t ) - x t I dt. (6) 



70 Wo / 



We study a nominal feasible trajectory to = (i, it). Next we present a qualifi- 
cation hypothesis that is assumed throughout the article. Consider the mapping 

G:R n xU -> R d " 

(x ,u) i-> J7(x ,xt), 

where is the solution of |2]) associated to (xq,u). 
Assumption 2.1. The derivative of G at (xq,u) is onto. 

Assumption |2 . 1 1 is usually known as qualification of equality constraints. 



Definition 2.2. It is said that the trajectory w is a weak minimum of prob- 
lem (P) if there exists e > such that w is a minimum in the set of feasible 
trajectories w — (x,u) G W satisfying 

\\x- xWoo < s, \\u — u\\oo < £• 
The following first order necessary condition holds for w. 
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Theorem 2.3. If w is a weak solution then there exists A = (/3,p) £ R dr "* x X* 
such that p is solution of the costate equation 

-p t = D x H[X](x t ,u t ,t), a.e. on [0,T], (8) 

with transversality conditions 

Po = -D Xo t[X](x ,x T ), (9) 
Pt = D Xt £[\}(xq,x t ), (10) 

and the stationarity condition 

D u H[X]{x u u t ,t)=0, a.e. on [0,T], (11) 

is verified. 

It follows easily that since the pre-Hamiltonian H is affine in all the control 



variables, (11) is equivalent to the minimum condition 



H[X](x t ,ut,t) = min H[X](x t ,v,t), a.e. on [0, T]. (12) 



In order words, the element (w, A) in Theorem 2.3 satisfies the qualified Pontrya 



gin Maximum Principle and A is a Pontryagin multiplier. It is known that the 



Assumption 2.1 implies the existence and uniqueness of multiplier. We denote 
this unique multiplier by A = ((3,p). 

Let the switching function <& : [0,T] — > M m '* be defined by 

$ t := D u H[X]{x u u u t) = {ptMx t ))Zi- (13) 



Observe that the stationarity condition (11) can be written as 



$ t = 0, a.e. on[0,T]. (14) 

3 Optimality System 

In this section we present an optimality system, i.e. a set of equations that 
are necessary for optimality. We obtain this system from the conditions in 
Theorem |2.3| above and assuming that the strengthened generalized Legendre- 
Clebsch condition (to be defined below) holds. 

Observe that, since H is affine in the control, the switching function $ 



introduced in ( 13 ) does not depend explicitly on u. Let an index i = 1, . . . , m, 
and (d Mi $/di A ^) be the lowest order derivative of $ in which Ui appears with a 
coefficient that is not identically zero on (0,T). In Kelley et al. [24! it is stated 
that Mi is even, assuming that the extremal is normal (as it is the case here 
since w satisfies the PMP in its qualified form). The integer := Mi/ 2 is 
called order of the singular arc. As we have just said, the control u cannot be 



retrieved from equation ( 11 ). In order to be able to express u in terms of (p, x) 
from equation 

<I> t = 0, a.e. on[0,T], (15) 
we make the following hypothesis. 
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Assumption 3.1. The strengthened generalized Legendre-Clebsch condition (see 
e.g. 

Kelley fESH and Goh JUJ/j holds, i.e. 

d 



du 



$t^0, on[0,T]. 



(16) 



Here, by X >~ we mean that the matrix X is positive definite. Notice that 
function $ is affine in u, and thus u can be written in terms of (p, x) from ( 15 ) 
by inverting the matrix in (16 1. Furthermore, due to the regularity hypothesis 
imposed on the data functions, u turns out to be a continuous function of time. 

Hence, the condition ( 15 ) is included in our optimality system and we can 
use it to compute u in view of Assumption |3.1| In order to guarantee the 



(17) 



stationarity condition ( 14 ) we consider the endpoint conditions 

$ T = 0, <i> = 0. 



Remark 3.2. We could choose another pair of endpoint conditions among the 
four possible ones: <f?o = 0> 3>t = 0, $o = and $x = 0, always including at 
least one of order zero. The choice we made will simplify the presentation of 
the result afterwards. 

Notation: Denote by (OS) the set of equations composed by Q-Q,(|8])- 
(pl, (pi, fT7|, i.e. the system 

m 

x t = y^^u i:t fi(x t ), a.e. on [0,T], 



(OS) 



r}j(x ,x T ) =0, for j = l } ...,d v , 
-fit = D x H[X](x t ,Ut,t), a.e. on [0,T], 
p a = -D Xo £[X](x ,x T ), 
Pt = D XT £[X\(x ,x T ), 
$t=0, a.e. on[0,T], 
<J>t = 0, $o = 0. 



Let us give explicit expressions for $ and $. Denote the Lie bracket of two 
smooth vector fields g, h: W 1 —> M. n by 



[5, h](x) := g'(x)h(x) - h'(x)g(x). 
Define A: W l+m M nXn (M.) and B: R n -)• X„ xm (M) by 



A(x, u) := ^2 Uifi(x), B{x)v := ^ vji(x), 



(18) 



(19) 



i=0 



i=i 



for every v £ K m . Notice that the ith. column of B(x) is fi(x). For (x,u) E W 
satisfying |2]), let Bi(x t ,u t ) G M nxm (M.) given by 

d 



B 1 (x t ,u t ) := A(x t ,u t )B(x t ) - —B(x t ). 



In view of (19) and ([20]), the expressions in (17) can be rewritten as 
$t=PtB(x t ), $>t = -PtBi(x u u t ). 



(20) 



(21) 
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4 Shooting Algorithm 



The aim of this section is to present an appropriated numerical scheme to solve 
system (OS). For this purpose define the shooting function 



S: D(S) := 



yn+d n , 



( T](x Q ,x T ) \ 
Po + D x J{\}(x ,x T ) 
Pt ~ D XT £[\](x ,x T ) 
p T B(x T ) 
p Bi(x ,u 



(22) 



/ 



where (x,u,p) isasolution of (I2]),(|8|,( 15 1 corresponding to the initial conditions 
(xq,pq), and with A := (/3,p). Here we denote either by (01,02) or 



an 



element of the product space A\ x A%. Notice that the control u retrieved from 



( 15 ) is continuous in time, as we have already pointed out after Assumption 



3.1 Hence, wc can refer to the value uq, as it is done in the right hand-side of 
(22). Observe that in a simpler framework having fixed initial state and no final 



constraints, the shooting function would depend only on p . In our case, since 
the initial state is not fixed and a multiplier associated with the initial-final 
constraints must be considered, S has more independent variables. Note that 
solving (OS) consists of finding v g D(6>) such that 



S(y) = 0. 



(23) 



Since the number of equations in (23) is greater than the number of unknowns, 



the Gauss-Newton method is a suitable approach to solve it. This algorithm 
will solve the equivalent least squares problem 



min >S(V)| 

i/£D(S) 1 W 1 



(24) 



At each iteration k, given the approximate values v k , it looks for A k that gives 
the minimum of the linear approximation of problem 



min \S(is k ) +S'(v k )A\ 

AeD(S) 1 



Afterwards it updates 



k+l 



(25) 



(26) 



In order to solve the linear approximation of problem (25) at each iteration k, 
we look for A fc in the kernel of the derivative of the objective function, i.e. A fe 
satisfying 

S'(v k ) T S'{v k )A k + SV) T S(*/-) = 0. (27) 

Hence, to compute direction A k the matrix S' (v k ) T S' (v k ) must be nonsingular. 
Thus, Gauss-Newton method will be applicable provided that S' (v) T S' (v) is 
invertible, where v :— (xq,Pq, $)■ Easily follows that S' (i>) T S' (v) is nonsingular 
if and only if S'(P) is one-to-one. Summarizing, the shooting algorithm wc 



propose here consists of solving the equation ( 23 ) by the Gauss- Newton method 



defined by (26)-(27|. 
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Since the right hand-side of system (23 1 is zero, the Gauss-Newton method 
converges locally quadratically if the function S has Lipschitz-continuous deriva- 
tive. The latter holds here given the regularity assumptions on the data func- 
tions. This convergence result is stated in the proposition below. See, e.g., 
Fletcher [T5] or Bonnans for a proof. 

Proposition 4.1. IfS'(i>) is one-to-one then the shooting algorithm is locally 
quadratically convergent. 

The main result of this article is to present a condition that guarantees the 
quadratic convergence of the shooting method near the optimal local solution 
(w, A). This condition involves the second variation studied in Dmitruk [20ll2~T] . 
more precisely, the sufficient optimality conditions therein presented. 

4.1 Linearization of a Differential Algebraic System 

For the aim of finding an expression of S'(u), we make use of the linearization 
of (OS) and thus we introduce the following concept. 

Definition 4.2 (Linearization of a Differential Algebraic System). Consider a 
system of differential algebraic equations (DAE) with endpoint conditions 

Ct=T(Ct,a t ) 1 (28) 
= G(tt,a t ), (29) 
= Z(Co,Ct), (30) 

where T : R m+n -> R n , Q : R m+n -> R d ° and 1 : R 2n -> R d ? are C 1 functions. 
Let (C , ot ) be a C 1 solution. We call linearized system at point ((°,a°) the 
following DAE in the variables £ and a, 

Ct = Lin J 7 | (£0,0,0) (Ct.at), (31) 

= Lin^ |( f o jQ .o) (Ct,a t ), (32) 

= LinZ l^o^oj (Co, Ct), (33) 

where 

Lin J 7 1(^0) (C t ,a t ) :=P((°,a° t )(( u a t ), (34) 
and the analogous definitions hold for LinC/ and Lin'H. 

The technical result below will simplify the computation of the linearization 
of (OS). Its proof is immediate. 

Lemma 4.3 (Commutation of linearization and differentiation). Given Q and 
J- as in the previous definition, it holds 

-^Lin£ = Lin ^-Q, — Lin J 7 = Lin ^-J 7 . (35) 
at at at at 
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4.2 Linearized optimality system 

In the sequel, whenever the argument of functions A, £?,£?!, etc. is omitted, 
assume that they are evaluated at the reference extremal (w,X). Define the 
m x n— matrix C, the n x n— matrix Q and the m x n— matrix M by 



C := H ux , Q := H xx , M := B Q -G-CA. 



(36) 



Notice that the zth. row of matrix C is the function pf[, for i = 1, ...,m. 
Denote with (z, v, A := (f3,q)) the linearized variable A = (/3,p)). In view 
of equations (21 ) and p6| we can write 



Lin 5> t =q t B t + zj Cj . 



(37) 



The linearization of system (OS) at point (£, m, A) consists of the linearized state 
equation 

z t = A t z t + B t v u a.e. on [0, T], 
with endpoint conditions 



= Dr)(x a ,x T )(zo,z T ), 
the linearized costate equation 

-q t = q t A t + zjQ t + vJC t , a.e. on [0, T] 
with endpoint conditions 



'7o 



q T 



zjDl 2 e + z^D 2 XoXT i + J2^ D -o 



Vj 



3 = 1 



z+DUl + z* Di oXT l- 



(x ,Xt) 



and the algebraic equations 



(38) 
(39) 

(40) 

(41) 
(42) 



= Lin $ = ——^{qB + Cz), a.e. on[0,T], 
= Lin $r = qxBx + CtZt, 

= Lin $ = -^:(q B + Cz )t=o- 
at 



Here we used equation (371 and commutation property of Lemma 4.3 
Q and p7| . Observe that (|43])-([47} and Lemma [4~3| yield 

= Lin $ t = q t B t + zJCj, on [0, T], 



(43) 
(44) 
(45) 
to write 

(46) 



and 



= Lin $ t = -qBi - z T M T + v 1 (-CB + B 1 C" ), a.e. on [0, T] 



T/->T\ 



By means of Theorem |5.2| to be stated in Section [5] afterwards we can see that 
the coefficient of v in previous expression vanishes, and hence, 



= Lin $ t = -qBi — z M , on [0, T] 



(47) 
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Note that both equations (46) and (47) hold everywhere on [0,T] since all the 
involved functions are continuous in time. 



Notation: denote by (LS) the set of equations (38 1-([47|). 

Once we have computed the linearized system (LS), we can write the deriva- 
tive of S in the direction D := (zo, <70j p) as follows. 



S'{v)v 



qa + 
q T - 



V 



Dr](x ,x T )(z , z T ) 

zlDl % i + zjD 2 X0XT £ + PjD XtVj 

q T B T + z^Cj 
q B lfi + z T M T 



(x ,x T ) 
(x ,x T ) 



(48) 



where (v, z, q) is the solution of |38|,(40|,(43) associated with the initial condi- 
tion (zo,<7o) and the multiplier j3. Thus, we get the property below. 



Proposition 4.4. S'{v) is one-to-one if the only solution of (38) -(40 1, (43) with 
the initial conditions zq = 0, qo — and with [3 = is (w, z, q) = 0. 



5 Second Order Optimality Conditions 

In this section we summarize a set of second order necessary and sufficient 
conditions. At the end of the section we state a sufficient condition for the local 
quadratic convergence of the shooting algorithm presented in Section [4j The 
latter is the main result of this article. 



Recall the matrices C and Q defined in ( 36 1 , and the space W given at the 



beginning of Section [2] Consider the quadratic mapping on W, 

r T 

Q(z,v) := \D 2 l{z ,z T f + \ \ [z T Qz + 2v T Cz] dt. (49) 

Jo 

It is a well-known result that for each (z, v) £ W, 

D 2 C{z,vf = fl(z,v). (50) 

We next recall the classical second order necessary condition for optimality that 
states that the second variation of the Lagrangian function is nonnegative on 
the critical cone. In our case, the critical cone is given by 



C := {(z,v) e W : (38)-d39| hold}, (51) 



and the second order optimality condition is as follows. 

Theorem 5.1 (Second order necessary optimality condition). If w is a weak 
minimum of (P) then 

Q(z, v) > 0, for all (z, v) € C. (52) 

A proof of previous theorem can be found in, e.g., Levitin, Milyutin and 
Osmolovskii [28]. The following necessary condition is due to Goh [26] and it 
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is a nontrivial consequence (not immediate) of Theorem 5.1 Define first the 
m x m— matrix 



R := B T QB - CB X - (CB 1 ) T 



df 



(CB). 



(53) 



Theorem 5.2 (Goh's Necessary Condition). If w is a weak minimum of (P), 
then 

CB is symmetric, (54) 

and 

RhO. (55) 



Remark 5.3. Observe that (54) is equivalent to pf-fj — pf'jfi, for every pair 
i,j = 1, . . . ,m. These identities can be written in terms of Lie brackets as 

P[fh fj] = °) for i,j = 1, • • ■ ,m. (56) 



Notice that (54 1 implies, in view of (53), that R is symmetric. The components 

p[fi,[fjJo]}, (57) 



of matrix R can be written as 



R, 



and hence, its symmetry implies 

p[fi,[fj,fo]]=p[fj,[fiJo]], for i,j = l,...,m. (58) 

The latter expressions involving Lie brackets can be often found in the literature. 

The result that we present next is due to Dmitruk [20 and is stated in terms 
of the coercivity of £1 in a transformed space of variables. Let us give the details 
of the involved transformation and the transformed second variation. Given 
(z,v) G W, define 



v s ds, 



£t ■= z t - B(x t )y t . 



(59) 
(60) 



This change of variables, first introduced by Goh in [25], can be perform in any 
linear system of differential equations, and it is known as Goh's transformation. 

We aim to perform Goh's transformation in ( [49] ). To this end consider the 
spaces U 2 := L 2 (0, T; R m ) and X 2 := W%(0, T; R n ), the function g : U 2n + m R 
defined by 



9(Co, Cr, h) := D 2 £ (Co, Ct + B T hf + h T C T {2( T + B T h), 



(61) 



and the quadratic mapping 
fi : X 2 x U 2 x W n -> K 



(t,y,h) ^y(to,fr,h) + i f {CQZ + 2y T MZ + y T Ry}dt, 

Jo 



(62) 



where the involved matrices where introduced in (19), (36) and (53). 
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Proposition 5.4. If w is a weak minimum of (P), then 

n(z,v) = &,(£,, y,y T ), 
whenever (z,v) G W and (£,y,y T ) G X x y x R m saiis/y (fggb-dgO >. 



(63) 



The latter result follows by integrating by parts the terms containing v in 



(49), and by replacing z by its expression in (60). See, e.g., Aronna et al. 
for the detailed calculations that lead to ( 63 ) . 

Define the order function 7:R"xM 2 x R m — > K as 



7(£o,2/,/i):=|£o| 2 + / y?dt+\h\ 



(64) 



We call (Sx, v) G W a feasible variation for w if (i + &c, u + t>) satisfies ([2])-(|3]). 

Definition 5.5. W^e say i/iai w satisfies the 7— growth condition in the weak 
sense if there exists p > such that, for every sequence of feasible variations 
{(5x k ,v k )} converging to in W, 



J{u + v k ) - J(u) > P1 ($,y k ,y k ), 
holds for big enough k, where y k := J Q w^'ds, and £ k is given by (60). 



(65) 



In previous definition, given that (Sx k ,v k ) is a feasible variation for each k, 
the sequence {(Sx k , v k )} goes to in W if and only if {v k } goes to in U. 



Observe that if (z, v) G W satisfies (|38|)-(39), then (£,y,h := y-r) given by 
transformation (59l-((60| verifies 



t = M + Biy, 

Dn{x ,x T )(S,o,t,T + B T h) =0. 



Set the transformed critical cone 

V 2 := {(£,y,h) EX 2 xU 2 x 



(66)-(67) hold}. 



(66) 
(67) 



(68) 



The following is an immediate consequence of the sufficient condition estab- 
lished in Dmitruk [20] (or [211 Theorem 3.1]). 

Theorem 5.6. The trajectory w is a weak minimum of (P) satisfying 7— growth 
condition in the weak sense if and only if (54) holds and there exists p > such 
that 

ft(t,y,h)>rY(ta,V,h), onP 2 . (69) 

The result presented in |20| applies to a more general case having finitely 
many equalities and inequalities constraints on the initial and final state, and a 
set of multipliers consisting possibly of more than one element. 



Remark 5.7. If (69) holds then necessarily 

R y pi, 



(70) 



where I represents the identity matrix. 
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Theorem 5.8. Ifw is a weak minimum of (P) satisfying (69), then the shooting 
algorithm is locally guadratically convergent. 

We present the proof of previous theorem at the end of Section [7| 



Remark 5.9. It is interesting to observe that condition (69) is a quite weak 



assumption in the sense that it is necessary for ~f— growth and its corresponding 



relaxed condition ( 52 ) holds necessarily for every weak minimum. 



Remark 5.10 (Verification of (69)). The sufficient condition in (69) can be 
sometimes checked analytically. On the other hand, when the initial point £o is 
fixed, it can be characterized by a Riccati-type equation and/or the nonexistence 
of a focal point as it was established in Zeidan 131]/ . Furthermore, under certain 



hypotheses, the condition ( 69 I can be verified numerically as proposed in by 
Bonnard, Caillau and Trelat (see also the survey in U3f ). 



6 Corresponding LQ Problem 

In this section we study the linear-quadratic problem (LQ) given by 

&(€,y,hT) min, 



([66J-(|67J, 

h = 0, ho free. 



(71) 
(72) 
(73) 



Here y is the control, £ and h are the state variables. Note that if condition 



(69) holds then (LQ) has a unique optimal solution (£, y, h) = 0. Furthermore, 
recall that (69) yields (70) as it was said in Remark 5.7 In other words, ( [69] ) 
implies that the strengthened Legendre-Clebsch condition holds at (£, y, h) = 0. 
Hence, the unique local optimal solution of (LQ) is characterized by the first 
optimality system, that we denote afterwards by (LQS). In Section[7]we present 
a one-to-one linear mapping that transforms each solution of (LS) (introduced 
in section 4.2) into a solution of this new optimality system (LQS). Theorem 



PI will follow. 

Denote by x an d Xh the costate variables corresponding to £ and h, respec- 
tively; and by (3 L Q the multiplier associated to the initial-final linearized state 
constraint (|67|). Note that the qualification hypothesis in Assumption 2.1 im- 



plies that {Dr)j(xo,XT)}j=i are linearly independent. Hence any weak solution 
(£, y, h) of (LQ) has a unique associated multiplier \ L ® := (x, Xh>P L ®) solution 
of the system that we describe next. The pre-Hamiltonian of (LQ) is 



n^ LQ K^y) ■= X(A£ + B x y) + i(£ T Q£ + 2y T Mi + y T Ry). 



(74) 



Observe that T-L does not depend on h since the latter has zero dynamics and 
does not appear in the running cost. The endpoint Lagrangian is given by 



i LQ [X LQ ](io^T,h T ) := y(^o^T,h T ) + J2^ Q Dv 3 ^o^t + B T h T ). (75) 

3=1 



The costate equation for x is 

-x = D^n[\ LQ ] = x a + CQ + y T M, 



(76) 
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with endpoint conditions 

Xo = -D io e L Q[x L Q] 



UDIJ + (Ct + B T h) 1 Dl oXT (. + J2%i Pf Q Dx Vj 



(77) 



xt= D iT e LQ [\ LQ ] 

= $Dl oXT l + (£ T + B T hYDl 2 £ + h^Cr + pf Q D XTVj . (78) 



For costate variable Xh we get the equation 

Xh = 0, 
Xhfi = o, 

Xh,T = D h l L Q[\ L % 



(79) 
(80) 
(81) 



Hence, \h = and thus (81) yields 



= $D 2 XoXT £B T + (fr + B T h) T {D 2 xl m T + Cf ) + £ ^ Q D XTm B T . (82) 

3=1 



The stationarity with respect to the new control y implies 



= DyH = x-Bi + C T M T + y T R 



(83) 



Notation: Denote by (LQS ) the set of equations consisting of ( 66 )- ( 67 ) , |73[) , ( 76 1 



(|78j) , (|82j> and (|83j), i.e. (LQS) is the system 

f £ = M + B lV , 

Dn{x Q ,x T )(£,o,tT + B T h) = 0, 

/i = 0, 



Xo = 



XT 



ijD 2 x U + (fr + B T h) T D 2 XoX J + ^ /^A^j 

3=1 

totiioxj + (fr + BThfB-l^l + h T C T + ]T pf Q D XT r, h 

3=1 

= $D 2 X0XT tB T + (fr + S t ^) T (^^Bt + C?) + X] P- Q D XTnj B T , 

(LQS) 



3 = 1 



I o = x s 1 + e T M T + ? / T i? 



Notice that (LQS) is a first order optimality system for problem ( 71 1-( 73 ) 



7 The Transformation 

In this section we show how to transform a solution of (LS) into a solution of 
(LQS) via a one-to-one linear mapping. Given (z, v , q, (3) e X x U x X* x M. dr "* , 
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define 



Vt 



v s ds, £ := z- By, x ■= q + V C, Xh ■= 0, h := y T , (3 



LQ 



Pi- (84) 



The next lemma shows that the point (£,y,h,x,Xh, (3 L ®) is solution of (LQS) 
provided that (z,v,q,/3) is solution of (LS). 



Lemma 7.1. The one-to-one linear mapping defined by (84) converts each so- 
lution of (LS) into a solution of (LQS). 



Proof Let (z,v,q,(3) be a solution of (LS), and set (S,,y,x, P LQ ) by 
Part I. We shall prove that (£, y, x, fi L ®) satisfies conditions ( [66] ) and (67). 
Equation (66) follows by differentiating expression of £ in ( p4| , and equation 
([67] follows from ([39]). 



Part II. We shall prove that ($;,y,X, P LQ ) verifies ((76])-((78j) and Differen- 
tiate x m (84), use equations ( |40| ) and ([84]), recall definition of M in (36) and 
obtain _ 

~X= -q-v T C -y T C 
= qA + z T Q - y T C 

= x A + eQ + y T (-CA + B T Q-C) 



(85) 



Hence (|76[) holds. Equations (77 1 and (78) follow from (41) and (42). Combine 




(|42J) and (|44|) to get 




q T B T + z± Cy. 



(i ,i T ) 



+ ZrpCj,. 



(86) 



Performing transformation ( 84 ) in the previous equation yields ( 82 ) . 
Part III. We shall prove that (83) holds. Differentiating (46) we get 

d 







d* 



-Lin $ 



dt 



(qB + z'C). 



(87) 



Consequently, by ( 38 ) and ( 40 ) , 



= -(qA + z l Q + v l C)B + qB + (z l A l + v 1 B 1 )C 1 + z T C T , 



where the coefficient of v vanishes in view of ( 54 1 . Recall ( 20 ) and ( 36 1 . Per- 
forming transformation ( 84 ) and regrouping the terms we get from ( 88 ) , 

= - X Bi - CM T + y T {CB 1 - B T QB + B T A T C T + B T C T ). (89) 



Equation (83) follows from |53[) and condition ([54]). 

Parts I, II and III show that (£, y, x, (3 L ® ) is a solution of (LQS), and hence 
the result follows. 

□ 

Remark 7.2. Observe that the unique assumption we needed in previous proof 
was Goh's condition (54) that follows from the weak optimality ofw. 
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Proof, [of Theorem 5.8 We shall prove that ([69]) implies that S'(P) is one-to- 



one. Take (z,v,q,f3) a solution of (LS), and let (£, y, x, Xhi fi L ®) be denned by 
(84), that we know by Lemma 7.1 is solution of (LQS). As it has been already 
pointed out at the beginning of Section [6] condition (69) implies that the unique 
solution of (LQS) is 0. Hence (£, y, x, Xh> (3 L ®) = and thus (z,v,q,/3) = 0. 
Conclude that the unique solution of (LS) is 0. The latter assertion implies, 
in view of Proposition 4.4 that S'(v) is one-to-one. The result follows from 
Proposition WA\ □ 



8 Control Constrained Case 

In this section we add the following bounds to the control variables 

0<«i, t <l, for a.a. fe [0,T], fori = l,...,m. (90) 
Denote with (CP) the problem given by Q-© and @. 



Definition 8.1. A feasible trajectory w € W is a Pontryagin minimum of ( CP) 
if for any positive N there exists en > such that w is a minimum in the set 
of feasible trajectories w — {x,u) £ W satisfying 

\\x - x\\oo < £jv> \\ u ~ u\\i < £ N, II" ~ ""Hoc < N. 

Given i = 1, . . . , m, we say that Ui has a bang arc in (a, b) C (0, T) if u^j = 
a.e. on (a, b) or — 1 a.e. on (a, 6), and it has a singular arc if < iii.t < 1 
a.e. on (a, b). 

Assumption 8.2. Each component Ui is a finite concatenation of bang and 
singular arcs. 

A time t € (0, T) is called switching time if there exists an index 1 < i < m 
such that Ui switches at time t from singular to bang, or vice versa, or from one 



bound in ( 90 ) to the other. 



Remark 8.3. Assumption \8 . S\ rules out the solutions having an infinite number 
of swit- chings in a bounded interval. This behavior is usually known as Fuller's 
phenomenon (see Fuller J331). Many examples can be encountered satisfying 
Assumption\8.£\ as is the case of the three problems presented in Section [J7j[ 



With the purpose of solving (CP) numerically we assume that the structure 
of the concatenation of bang and singular arcs of the optimal solution w and 
an approximation of its switching times are known. This initial guess can be 
obtained, for instance, by solving the nonlinear problem resulting from the dis- 
cretization of the optimality conditions or by a continuation method. See Bctts 
[3"4"j or Biegler [3S] for a detailed survey and description of numerical methods 
for nonlinear programming problems. For the continuation method the reader 
is referred to Martinon [9]. 

This section is organized as follows. From (CP) and the known structure 
of u and its switching times we create a new problem that we denote by (TP). 
Afterwards we prove that we can transform w into a weak solution W of (TP). 



Finally we conclude that if W satisfies the coercivity condition (69), then the 



shooting method for problem (TP) converges locally quadratically. In practice, 
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the procedure will be as follows: obtain somehow the structure of the optimal 
solution of (CP), create problem (TP), solve (TP) numerically obtaining W, 
and finally transform W to find w. 

Next we present the transformed problem. 

Assumption 8.4. Assume that each time a control Ui switches from bang to 
singular or vice versa, there is a discontinuity of first kind. 

Here, by discontinuity of first kind we mean that each component of u has 
a finite nonzero jump at the switching times, and the left and right limits exist. 

By Assumption |8 . 2| the set of switching times is finite. Consider the partition 
of [0, T] induced by the switching times: 

{0 =: f < fx < . . . < fW_ x < f N := T}. (91) 

Set J fe := [f fc _i, f k ], and define for k = 1, . . . , N, 

Sk '■= {1 < i < m '■ ^ is singular on Ik}, (92) 
Eh := {1 < i < m : in — a.e. on Ik}, (93) 
Nk ■= {1 < i < m : Ui — 1 a.e. on Ik}- (94) 

Clearly S k U E k U N k = {1, . . . , m}. 

Assumption 8.5. For each k = 1, . . . , N, denote by us k the vector with compo- 
nents Ui with i 6 Sk- Assume that the strengthened generalized legendre-Clebsch 
condition holds on Ik, i.e. 

-^—H USh yO, on I k . (95) 
ou Sk k 

Hence, u$ k can be retrieved from equation 

H USk = 0, (96) 

since the latter is affine on us k as it has been already pointed out in Section [3] 
Observe that the expression obtained from ( 96 ) involves only the state variable 
x and the corresponding adjoint state p. Hence, it results that u$ k is continuous 
on Ik with finite limits at the endpoints of this interval. As the components Ui 
with i ^ Sk are either identically 1 or 0, we conclude that 

u is continuous on Ik- (97) 



By Assumption 8.4 and condition (97) (derived from Assumption |8.5[ ) we 
get that there exists p > such that 

p < u,i,t < 1 — p, a.e. on Ik, for k = 1, . . . , N, i € Sk- (98) 

Next we present a new control problem obtained in the following way. For each 
k = 1, . . . , N, we perform the change of time variable that converts the interval 
Ik into [0, 1], afterwards we fix the bang control variables to their bounds and 
finally, we associate a free control variable to each index in Sk- More precisely, 
consider for k — 1, . . . , N the control variables u\ g Too(0, 1; K), with i e Sk, 
and the state variables x k € H^(0, 1; E"). Let the constants Tk G K, for k 
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1, . . . , N — 1, which will be considered as state variables of zero-dynamics. Set 
Tq := 0, T/v := T and define the problem on the interval [0, 1] 



<p Q (xl,Xi) -> min, 

x k = (T fe - Tjk-i) I 2 

\iew fc u{o} 

T fe = 0, fc = l,...,JV-l, 



fc = 1,. 



,JV-1. 



(99) 

.,N, (100) 

(101) 
(102) 
(103) 



Denote by (TP) the problem consisting of equations (99)-( T03| . The link be- 
tween the original problem (CP) and the transformed one (TP) is given in 
Lemma 8.6 below. Set for each k — 1, . . . , N : 



:= x{f u _ x + (ft - f fe _ 1 )s), for s E [0, 1], 

:= Uj(T fe _i + (T k - T fe _i)s), for i e S k , a.a. s € [0, 1]. 



Set 



fc\iV 



*JV-1\ 



(104) 
(105) 

(106) 



Lemma 8.6. is a Pontryagin minimum of ( CP), then W is a weak solution 
of (TP). 

Proof. The idea of the proof is to derive the weak optimality of W from the 
Pontryagin optimality of w and condition ( 98 ) . Since w is a Pontryagin mini- 
mum for (CP), there exists e > such that w is a minimum in the set of feasible 
trajectories w = (x,u) satisfying 

\\x - x\\oo < £, — u||l<£, ||lt- filloo < 1. (107) 

Consider 8, e > 0, and a feasible solution ((x k ), (u k ), (Tk)) for (TP) such that 
\T k -f k \<8, \\u k - ufWn < e, for all k = 1, . . . , N. (108) 

, N. Denote 



We shall relate s in ( |107| with 8 and e in ( 108| ). Let k = 1, 
I k := (Tfe_x, Tfe), and define for each i = 1, . . . , m : 





,fe ( t-T fc _x \ 
1, 



if t e I k and z € 
if £ € Ife and i € S k , 
i{ t e h and i e N k . 



(109) 



Let # be the solution of |2]) associated to it and having x = xj. We shall prove 
that (a;, it) is feasible for the original problem (CP). Observe that condition 

( |103| ) implies that x t — x k J when t £ I k , and thus £i = a;^. It follows 

that ([3]) holds. We shall check condition ([90]). For i e E k U N k , it follows from 
the definition in (109 1. Consider now i £ S k . Since (98) holds, by (105) we get 



p < u k s < 1 — p, a.e. on (0, 1). 



(110) 
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Thus, by ( 108 1 and if e < p, we get < u\ s < 1 a.e. on (0, 1). This yields 

< u ht < 1, a.e. on I k , (111) 

and thus the feasibility of (x,u) for (CP). 

We now estimate llu — u\\i. For k = 1, . . . ,N and i E Sk, 



dt 



t—Tk-i 
Tk — Tk-i 



dt. 



(112) 



Note that by Assumption 8.4 and condition (97), each uf is uniformly continuous 



on Ik, and thus there exists 9ki > such that if |s— s'\ < 9ki then s — u^ s , \ < 
s. Set 9 :— min Of-i > 0. Consider then 5 such that if \Tk — Tk\ < 5, then 

t—Tk-i t—Tk-i 



< 0. From (1081 and (112 1 we get 



\ui it - Ui,t\dt < 2emeas(/ fc n 4). 



Iknl k 



Assume, w.l.o.g., that Tk < Tfc and note that 

,. I t-Tk-i 



\ui, t -Ui,t\dt < 



Tk — Tk-i 



t — Tk-\ 
Tk — Tk-\ 



(113) 



dt<5e, (114) 



where we used (108) in the last inequality. From (113) and (114) we get \\ui 
Ui\\x < e(2T+ (N- 1)6). Thus < e if 



e(2T+ (N — 1)5) < e/m. 



(115) 



We conclude from ( 107 ) that ((x k ), (uf), (Tk)) is a minimum on the set of feasible 
points satisfying (108) and ( |115 ). Thus FT is a weak solution of (TP), as it was 
to be proved. □ 

We shall next propose a shooting function associated to (TP). The pre- 
Hamilto- nian of the latter is 



H := J2( T k ~ T k -i)H k 



(116) 



k=i 



where, denoting by p k the costate variable associated to x k , 



\iGAf fc U{0} iSSfc 

Observe that Assumption |8.5| made on u yields 



~®-h u ^0, on [0,1], 
ou 



(117) 



(118) 



i.e. the strengthened generalized Legendre-Clebsch condition holds in problem 
(TP) at w. Hence we can define the shooting function for (TP) as it was done 
m Section |4] for (P). 
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The endpoint Lagrangian is 



JV-l 



l:= +T,PM( x 0' x ^ + E ^ 



„k+ls 



3 = 1 



k=l 



The costate equation for p k is given by 



— {Tk — T k _i)D x kH k , 



with endpoint conditions 



-D„i 



Pi 
Po 



i=i 

for fc = 1, ... ,7V - 1, 
for k = 2, . . . , N, 

i = D x Nip + PjD x xrr)j. 

3 = 1 



1\ 

q/s-1 



For the costate variables p Tk associated with we get the equations 



-H" + H 



k+l 



0, pf fc =0, for k = l,...,N- 1, 



(119) 
(120) 

(121) 
(122) 
(123) 

(124) 



Remark 8.7. We can sum up the conditions in (124 1 integrating the first one 
and obtaining J Q (H k+1 — H k )dt = 0, and hence, since H k is constant on the 
optimal trajectory, we get the equivalent condition 



H k = H { 



k+l 



for k = 1, 



,JV-1. 



(125) 



So we can remove the shooting variable p Tk and keep the continuity condition 
on the pre-Hamiltonian. 



Observe that ( 103 ) and ( 122 ) imply the continuity of the two functions ob- 
tained by concatenating the states and the costates, i.e. the continuity of X 
and P defined by 

X :=xl, X s :=x k {s-{k-l)) 1 for s € (k- l,k], k= 1,...,N, (126) 
P :=pl, P s :=p k (s-(k-l)), forse(fc-l,fc], k = l,...,N. (127) 

Thus, while iterating the shooting method, we can either include the condi- 
tions ( 103 ) and ( 122 ) in the definition of the shooting function or integrate the 
differential equations for x k and p k from the values x\~ l and p\ _1 previously 
obtained. The latter option reduces the number of variables and hence the size 
of the problem, but is less stable. We shall present below the shooting function 
for the more stable case. For this end define the n x n~ matrix 



ieN k u{o} ies k 
the n x |Sfe|— matrix B k with columns fi{x k ) with i £ Sk, and 



(128) 



B? := A k B k - —B 
1 At 



(129) 
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We shall denote by gi(x k ,u k ) the ith. column of B k for each i in Sk- Here u k is 
the | Sk\— dimensional vector of components u k . The resulting shooting function 
for (TP) is given by 



S: 



vNn+N-l 



|d„+(JV-l)n x 
/ 



((x k ),(T k ),(p k ),[3)=:v h>S(i/) 



(4 - 4 + Vi,..,iV-l 
pl + D x {i[\]{xl.x^) 

(p k -p% +1 ) k =l,...,N-l 

p?-D x «i[\](xlx?) 

(,-"1 — -"0 )k=l,...,N-l 
(Po/i(4))fe=l,...,JV, *6S fc 

V (Po&^o^o))*^!,...,^, *es fc 



(130) 

Here we put both conditions H u = and H u — at the beginning of the interval 
since we have already pointed out in Remark |3.2| that all the possible choices 
were equivalent. 

Since problem (TP) has the same structure than problem (P) in section [2] 
i.e. they both have free control variable (initial-final constraints), we can apply 
Theorem |5.8| and obtain the analogous result below. 

that w is a Pontryagin minimum of ( CP) such that W 



Theorem 8.8. At 



defined in (106) satisfies condition (69) for problem (TP). Then the shooting 
algorithm for (TP) is locally quadratically convergent. 



Remark 8.9. Once system (130) is obtained, observe that two numerical im- 



plementations can be done: one integrating each variable on the interval [0, 1] 
and the other one, going back to the original interval [0,T], and using implicitly 
the continuity conditions (103), (122) and (125) at each switching time. The 
latter implementation is done in the numerical tests of Section \10\ below. In 
this case, the sensibility with respect to the switching times is obtained from the 
derivative of the shooting function. 



8.1 Reduced Systems 

In some cases we can show that some of the conditions imposed to the shooting 



function in (130) arc redundant. Hence, they can be removed from the formu- 



lation yielding a smaller system that we will refer as reduced system and which 
is associated to a reduced shooting function. 

Recall that when defining S we are implicitly imposing that H u = 0. The 
latter condition together with H u = H u i = 0, both included in the definition 
of S, imply that H u = H u = 0. Hence, 



p\h{x k ) = p k g t (x k ,u k ) =0, for k = 1, . . . , N, i e S k 



(131) 



and, in view of the continuity conditions (103) and (122) 

l M4 +1 ) 



P k +1 



p k+1 gi (x k+1 ,u k+1 ) = 0, forfc=l,. 



N —1, ieS k . (132) 
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Therefore, if a component of the control is singular on I k and remains being 
singular on I k +i, then there is no need to impose the boundary conditions on 

H u and H u since they are a consequence of the continuity conditions and the 
implicit equation H u = 0. 



( 132 1 we obtain, 



Observe now that from (117), (130) and previous two equations (131) and 



N k U{0} N k u{Q}\S k+1 

On the other hand, 

<' +1 =Po fc+1 E /^o +1 )- (134) 

JV fc + 1 u{0}\S fc 

Thus, fff = #o +1 if N k U {0}\Sfc + i = N k+1 U {0}\S k . The latter equality 
holds if and only if at instant T k all the switchings are either bang-to-singular 
or singular-to-bang. 

Definition 8.10 (Reduced shooting function). We call reduced shooting func- 



tion and we denote it by S r the function obtained from S defined in ( 130 1 by 
removing the condition H\ = Hq +1 whenever all the switchings occurring at T k 
are either bang-to- singular or singular-to-bang, and removing 

pg/i(a:g) = 0, p§fli(asg,«g)=0, (135) 

for k = 2, . . . ,N and i G S k -i H S k - 

8.2 Square Systems 

The reduced system above-presented can occasionally result square, in the sense 
that the reduced function S r has as many variables as outputs. This situation 



occurs, e.g., in problems 1 and 3 of Section [10] The fact that the reduced 
system turns out to be square is a consequence of the structure of the optimal 
solution. In general, the optimal solution u yields a square reduced system if 
and only if each singular arc is in the interior of [0, T] and at each switching time 
only one control component switches. This can be interpreted as follows: each 
singular arc contributes to the formulation with two inputs that are its entry and 
exit times, and with two outputs that correspond to Po/j( x o) = 9i{ x o> u o) = 0' 
being I k the first interval where the component is singular and i the index of the 
analyzed component. On the other hand, whenever a bang-to-bang transition 
occurs, it contributes to the formulation with one input for the switching time 
and one output associated to the continuity of the pre-Hamiltonian (which is 
sometimes expressed as a zero of the switching function). 

9 Stability under Data Perturbation 

In this section we investigate the stability of the optimal solution under data 



perturbation. We shall prove that, under condition (69), the solution is stable 



under small perturbations of the data functions po, fi and 77. Assume for this 
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stability analysis that the shooting system of the studied problem can be reduced 
to a square one. We gave a description of this situation in Subsection |8.2| Even 
if the above-mentioned square systems appear in control constrained problems, 
we start this section by establishing a stability result of the optimal solution for 
an unconstrained problem. Afterwards, in Subsection |9.2[ we apply the latter 
result to problem (TP) and this way we obtain a stability result for the control 
constrained problem (CP). 



9.1 Unconstrained control case 



Consider then problem (P) presented in Section [2j and the family of problems 
depending on the real parameter fi given by: 



iPo(x , x T ) -> min, 

m 

x t = J2 u i,tft( x t)i forte (0,T), 

t=0 

77 /i (x ,x T ) = 0. 



(Pm) 



Assume that <p% : K 2n+1 — > K and rf : E 2 ™ +1 — > R d " have Lipschitz-continuous 
second derivatives in the variable (xq,xt) and continuously diffcrcntiable with 
respect to /i, and /f : E™ +1 — > R" is twice continuously diffcrcntiable with 
respect to x and continuously differentiable with respect to the parameter fi. In 
this formulation, the problem (Po) associated to /i = coincides with (P) , i.e. 
¥o — ^Oi fi — fi f° r * — 0, . . . , m and rf = r\. Recall (69 1 in Theorem 5.6 and 
write the analogous condition for (P M ) as follows: 



W(li,y,h)>fr/(t; ,y,h), 



on Vg, 



(136) 



where and Vi? are the second variation and critical cone associated to (P A 
respectively. Let <S M be the shooting function for (P^). Thus, we can write 



(137) 



where we indicate with M the dimension of the domain of S. The following 
stability result will be established. 

Theorem 9.1 (Stability of the optimal solution). Assume that the shooting 
system generated by problem (P) is square and let w be a solution satisfying the 
uniform positivity condition ( 69 1 . Then there exists a neighborhood J C K of 



0, and a continuous differentiable mapping fi M> = (x^jU^), from J to W, 
where is a weak solution for ( P^ ) . Furthermore, verifies the uniform 
positivity (136). Therefore, in view of Theorems 5.6 and 5.8. the 7— growth 



holds, and the shooting algorithm for (P 11 ) is locally quadratically convergent. 

Let us start showing the following stability result for the family of shooting 
functions 



Lemma 9.2. Under the hypotheses of Theorem 9.1, there exists a neighborhood 
T C ffi of and a continuous differentiable mapping fi t— > - - (x^ ,Pq , /3^), 
froml to R M , such that S^{v^) — 0. Furthermore, the solutions {x^ ,p^) of 
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([2])-([8])-( 15 1 with initial condition (xq,Pq) and associated multiplier fj^ provide 
a family of feasible trajectories u> M := (x M , u^ 1 ) verifying 



i\\oo + IK - &IL + \\p» -pIIoo + 1/^ - /3| = o(m)- 



(138) 



Proof. Since (69) holds, the result in Theorem |5 .8| yields the non-singularity of 
the square matrix D^S°(v). Hence, the Implicit Function Theorem is applicable 
and we can then guarantee the existence of a neighborhood B C M M of v, a 
neighborhood I c K of 0, and a continuously differentiable function T :X B 
such that 

S^(T(p)) = 0, for all fi el. (139) 

Finally, write v* 1 :— T(/x) and use the continuity of DT on 1 to get the first part 
of the statement. 

The feasibility of w M holds since equation ( 139 1 is verified. Finally, the 
estimation ( 138 1 follows from the stability of the system of differential equation 
provided by the shooting method. □ 

Once we obtained the existence of this w M feasible for (P M ), we may wonder 
whether it is locally optimal. For this aim, we shall investigate the stability of 
the sufficient condition ( [69] ). Denote by f2 M and V% the quadratic mapping and 
critical cone related to (P M ), respectively. Given that all the functions involved 
in fV' are continuously differentiable with respect to fi, the mapping Cl^ itself 
is continuously differentiable with respect to (i. For the perturbed cone we get 
the following approximation result. 



Lemma 9.3. Assume the same hypotheses as in Theorem 9.1 Take /i € 1 and 
(^,y tl ,h fJ ') e V$. Then there exists (£,y,h) € V 2 such that 



£ \ + \\y> M -yh + \h»-h\ = O(L l ). 



(140) 



The definition below will be useful in the proof of previous Lemma. 
Definition 9.4. Define the function fj :U x W L — > M. " , given by 

fj(u,x ) := n(x ,x T ), 
where x is the solution of ^ associated to (it, xq). 



(141) 



Proof, [of Lemma [9 .3| Recall that Drj(u,xo) is onto by Assumption 2.1 Call 



back the definition of the critical cone C given in (511, and notice that we 
can rewrite it as C = {(z, v) € W : Q{z,v) = 0} = KerC*, with Q(z,v) :— 
Dr)(x , xt){zq, zt) being an onto linear application from W to R d, K In view of 
Goh's Transformation d59|> - d60b , 



(142) 



Di](x ,x T )(z ,z T ) = Dt](x ,x t )(£o,£t + B T y T ), 



for (z. v) £ W and (£, y) being its corresponding transformed direction. Thus, 
the cone Vi can be written as V2 = {C G H '■ ^-(C) = 0} = Ker/C, with 
C := (£,»,/»), H := X 2 x U 2 x R n , and /C(C) := Drj(x , z T )(£o, 6r + B T h). Then 
JC £ £(H,R d ") and it is surjective. Analogously, V% = {QeH: JC^C) = 0} = 
Ker/C", with 



(143) 
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Let us now prove the desired stability property. Take £ M E V2 = Ker /C M having 
- 1. Hence /C(C M ) = /C^C) + QC - /C^)(C M ), and by estimation ( fl43| ), 

TOI = OW. (144) 

Ker/C © Im/C T , there exists C M '* G ?T such that 

C := C +IC T (C'*) G Ker/C. (145) 



ICIIw 



Observe that, since % 



This yields = /C(C) =/C(C M ) + /C/C T (C Al '*) = (£ -JC^){^) + /C/C T (C M '*). Given 
that /C is onto, the operator JCKJ is invertible and thus 

= -(/C/C T )- 1 (/C-/C' t )(C' i ). (146) 

The estimation ( |144| ) above implies HC^IIw* = C(m)- It follows then from ( |145[ ) 
that — (\\-h = O(fi), and therefore, the desired result holds. □ 



Proof, [of Theorem |9.1| We shall begin by observing that Lemma 9.2 provides 
a neighborhood 1 and a class of solutions {(x^, u^ 1 ,^, (3^)}^! satisfying ( 138 1. 
We shall prove that — {x^\u^) satisfies the sufficient condition (136) close 
to 0. 

Suppose on the contrary that there exists a sequence of parameters — > 
and critical directions , y* 1 " , W k ) G T 5 ^ with 7(^ fc , y^ k , /i^ fe ) = 1, such that 



Since f2^ is Lipschitz-continuous in from previous inequality we get 

Cl(^ k ,y IJ "',h flk ) < o(l). 



In view of Lemma 9.3 there exists for each k, a direction (£ k ,y k 7 h k ) 
satisfying 

I Co - £0* I + 11/ - V^h + \h k - h>*\ = 0(ji k ). 



(147) 

(148) 
G V 2 

(149) 
(150) 



Hence, by inequality ([148J and given that w satisfies (691, 

Pi{£,ly\h k )<n(e,y\h k )<o{i). 

However, the left hand-side of this last inequality cannot go to since (£q , y k , h k ) 
is close to (^Q fc , y^ k , h^ k ) by estimation (149), and the elements of the latter 
sequence have unit norm. This leads to a contradiction. Hence, the result 
follows. □ 



9.2 Control constrained case 

In this paragraph we aim to investigate the stability of the shooting algorithm 
applied to the problem with control bounds (CP) studied in Section [8j Observe 
that previous Theorem |9.1| guarantees the weak optimality for the perturbed 
problem when the control constraints are absent. In case we have control con- 
straints, this stability result is applied to the transformed problem (TP) (given 



by equations ( 99 )-( 103 ) of Section |8j) yielding a similar stability property, but 
for which the nominal point and the perturbed ones are weak optimal for (TP). 
This means that they are optimal in the class of extremals having the same con- 
trol structure, and switching times and singular arcs sufficiently close in Lrx,. 
An extremal satisfying optimality in this sense will be called weak- structural 
optimal, and a formal definition would be as follows. 
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Definition 9.5 (Weak-structural optimality) . A feasible solution w for problem 
(CP) is called a weak- structural solution if its transformed extremal W given by 
(p4|-([l06j) is a weak solution of (TP). 



Theorem 9.6 (Sufficient condition for the extended weak minimum in the 
control constrained case). Let w be a feasible solution for (CP) satisfying As- 
sumptions 8.2 and 8.4 Consider the transformed problem (TP) and the corre- 



sponding transformed solution W given by (104) -(106 1. If w satisfies (69) for 



(TP), then w is an extended weak solution for (CP). 



Proof. It follows from the sufficient condition in Theorem 5.6 applied to (TP) 



□ 



Consider the family of perturbed problems given by: 
<p£(xo,x T ) -> min, 



it = £) «<,*//*(**)> forie(0,T), 

i=0 

n»(x ,x T ) = 0, 
0<ut<l, a.eon(0,T). 



(CP M ) 



The following stability result follows from Theorem 9.1 



Theorem 9.7 (Stability in the control constrained case). Assume that the re- 
duced shooting system generated by problem ( CP) is square. Let w be the solution 
of (CP) and {Tk}^ = i its switching times. Denote by W its transformation via 



equation (106). Suppose that W satisfies uniform positivity condition |69|) for 
problem (TP). Then there exists a neighborhood J C K of such that for every 
parameter fi G J there exists a weak- structural optimal extremal of (CP' 1 ) 
with switching times {Tt}^—^ satisfying the estimation 



N 

Ei T * 

k=l 



where It := 



T, 



N 

EE 

k=l ieSfc 



*»ii<x>,/?nJ* 



51100=00*), (151) 



(Tfc_i,Tk). Furthermore, the transformed perturbed solution W 1 



verifies uniform positivity ( 136 ) and hence quadratic growth in the weak sense for 
problem (TP) holds, and the shooting algorithm for (CP^ ) is locally quadratically 
convergent. 



9.3 Additional analysis for the scalar control case 

Consider a particular case where the control u is scalar. The lemma below shows 



that the perturbed solutions are Pontryagin extremals for ( CP^ ) provided that 
the following assumption holds. 

Assumption 9.8. (a) The switching function H u is never zero in the interior 
of a bang arc. Hence if u — 1 on (£1,^2) then H u < on (^1,^2)1 an d ifu = ~ 1 
on (tijta) then H u > on (ti,t2)- 

(b) If Tk is a bang-to-bang switching time then H u (Tk) =/= 0. 



The property (a) is called strict complementarity for the control constraint. 



RR n° 7763 



28 



Aronna & Bonnans & Martinon 



Lemma 9.9. Suppose that u satisfies Assumptio n\9.£^ Let 
9. 1 above. Then is a Pontryagin extremal for ( CP M I . 



;^ as in Theorem 



Proof. We intend to prove that w M satisfies the minimum condition ( 12 ) given 
by the Pontryagin Maximum Principle. Observe that on the singular arcs, 
H£ = since w M is the solution associated to a zero of the shooting function. 
It suffices then to study the stability of the sign of H£ on the bang arcs around 
a switching time. First suppose that u has a bang-to-singular switching at 
T k . Assume, without loss of generality, that u e 1 on 4 and u is singular on 
[T k ,T k+ i]. Let us write 

ff£ = a" + ti"^, (152) 

where a" and V* := -§^H£ are continuous functions on [0,T], and continuously 
diff erentiable with respect to fi since they depend on x^ and p^. Assumption 



5.5 



yields b° < on [T k ,T k+1 ], and therefore 

b»<0, on {T^T^}. 



(153) 



Due to (1521, the sign of H£ around depends on u^{Tj!+) - u»(Tg-). But 



this quantity is negative since u M passes from its upper bound to a singular arc. 



From the latter assertion and (1531 follows 



■) < 0, (154) 



Since H£ is null on [T£,T£ +1 ] 



and thus H£ is concave at the junction time T£. 
its concavity implies that it has to be negative before entering this arc. Hence, 
w 1 * respects the minimum condition on the interval I k . 

Consider now the case when u has a bang-to-bang switching at T k . Let us 
begin by showing that H£(T£) = 0. Suppose on the contrary that H£(T£) ^ 0. 
Then H^(T£+) — H^lTj?—) ^ 0, contradicting the continuity condition imposed 
on H in the shooting system. Hence H£{T?) = 0. On the other hand, since 
H u {T k ) ^ by Assumption 9.8 the value H£(T£) has the same sign for small 
ix. This implies that H£ has the same sign before and after that H u (before 
and after Tk), respectively. The result follows. □ 

Remark 9.10. We end this analysis by mentioning that if the transformed 
solution W satisfies the uniform positivity (69) for (TP), then w verifies the 
sufficient condition established in Aronna et al. JESj and hence it is actually a 



Pontryagin minimum. This follows from the fact that in condition ( 69 1 we are 
allowed to perturb the switching times, and hence ( 69 ) is more restrictive ( or 
demanding) that the condition in J 3 Of . 



10 Numerical Simulations 

Now we aim to check numerically the extended shooting method described 
above. More precisely, we want to compare the classical n x n shooting for- 
mulation to an extended formulation with the additional conditions on the pre- 
Hamiltonian continuity. We test three problems with singular arcs: a fishing 
and a regulator problem and the well-known Goddard problem, which we have 
already studied in [3S1 [37] . For each problem, we perform a batch of shoot- 
ings on a large grid around the solution. We then check the convergence and 
the solution found, as well as the singular values and condition number of the 
Jacobian matrix of the shooting function. 
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10.1 Test problems 
10.1.1 Fishing problem 

The first example we consider is a fishing problem described in |38) . The state 
xt G K represents the fish population (halibut), the control Ut € K is the fishing 
activity, and the objective is to maximize the net revenue of fishing over a fixed 
time interval. The coefficient (E — c/x) takes into account the greater fishing 
cost for a low fish population. The problem is 

T 

(E - c/x t ) u t U m ^dt, 

X t = TX t (1 - X t /k) - U t C/max, (Pi) 

0< Mt <l, Vt€[0,T], 
x = 70, xt free, 

with T=10,E = l,c= 17.5, r = 0.71, k = 80.5 and U max = 20. 

Remark 10.1. The state and control were rescaled by a factor 10 6 compared 
to the original data for a better numerical behavior. 

Remark 10.2. Since we have an integral cost, we add a state variable to adapt 
(Pi) to the initial-final cost formulation. It is well-known that its corresponding 
costate variable is constantly equal to 1. 

The pre-Hamiltonian for this problem is 

H := (c/x - E) u U max + p[r x (1 - x/k) -uU mllx \, (155) 
and hence the switching function 

$ t = D u H t = U max (c/x t -E-pt), Vte [0, T]. (156) 
The optimal control follows the bang-bang law 



u* t =0 if $ t > 0, 
u* t =1 if $ t < 0. 



(157) 



Over a singular arc where $ = 0, we assume that the relation $ = gives the 
expression of the singular control (t is omitted for clarity) 



t k r fee 2px 2px 2 

" singular = 2(c/x-p)U max [x~k~ P+ ^ W 



(158) 



The solution obtained for (Pi) has the structure bang-singular-bang, as 
shown on Figure [T] 

Shooting formulations. Assuming the control structure, the shooting un- 
knowns are the initial costate and the limits of the singular arc, 

v ■■= (po,ti,t 2 ) eM 3 . 

The classical shooting formulation uses the entry conditions on ^ 

Si(u) := (pT,$ tl ,$ tl ). 
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Figure 1: Fishing Problem 



Solving S\(y) — is a square nonlinear system, for which a quasi- Newton 
method can be used. Note that even if there is no explicit condition on t2 in S, 
the value of px does depend on ti via the control switch. 

The extended shooting formulation adds two conditions corresponding to the 
continuity of the pre-Hamiltonian at the junctions between bang and singular 
arcs. We denote [H] t := H t + — H t - the pre-Hamiltonian jump, and define 

&(!/) = (p 10 ,$ tl ,$ tl ,[H} tl ,[H} t2 ). (159) 

To solve Si(v) = we use a nonlinear least-square algorithm (see paragraph 



10.2 below for more details). 



10.1.2 Regulator problem 

The second example is the quadratic regulator problem described in Aly [3"9"] . 
We want to minimize the integral of the sum of the squares of the position and 
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speed of a mobile over a fixed time interval, the control being the acceleration. 



o 



x l,t — %2,ti 
%2,t = Ut, 

-l<u t <l, a.e. on[0,T] 
xq = (0, 1), xt free, 
T = 5. 



The corresponding prc-Hamiltonian 

H := ~{x\ + xf) +p\x 2 +P2U, 
and hence we have the switching function 

$ t := D u H t = p 2 ,f 
The bang-bang optimal control satisfies 

u* t = -sign p 2 , t if *t ^ 0. 
The singular control is again obtained from $ = and verifies 

^singular,* — x l,t- 



(Pa) 



(160) 



(161) 



(162) 



(163) 



The solution for this problem has the structure bang-singular, as shown on 
Figure [2] 

Shooting formulations. Assuming the control structure, the shooting un- 
knowns are 

v~ (Pi,o,J*,o,*i) eM3 - ( 164 ) 
For the classical shooting formulation, in order to have a square system, we can 
for instance combine the two entry conditions on $ and $, since we only have 
one additional unknown which is the entry time t\. Thus we define 

Siiv) := (pi,T,P2,T,$ 2 tl (165) 
The extended formulation does not require such a trick, we simply have 

S 2 {v) ■= (pi,T,P2,T,®ti,K,[ H ]ti)- (166) 



10.1.3 Goddard problem 

The third example is the well-known Goddard problem, introduced in Goddard 
pf0] and studied for instance in Seywald-Cliff [41] . This problem models the 
ascent of a rocket through the atmosphere, and we restrict here ourselves to 
vertical (unidimensional) trajectories. The state variables are the altitude, speed 
and mass of the rocket during the flight, for a total dimension of 3. The rocket 
is subject to gravity, thrust and drag forces. The final time is free, and the 
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STATE CONTROL 





Figure 2: Regulator Problem 



objective is to reach a certain altitude with a minimal fuel consumption, i.e. a 
maximal final mass. 



max m,T, 
r = v, 

v = -1/r 2 + l/m(T max? i - D(r,v)) 

fa = -&T max u, 

< Ut < 1, a.c. on (0, 1), 

r = 1, v = 0, to = 1, 

n = i.oi, 

T free, 



with the parameters b — 7, T max = 3.5 and the drag given by 

D(r,v) ^SlO^e" 500 ^- 1 ). 
The pre-Hamiltonian function here is 

H :=p r v+p v (-l/r 2 + l/TO(T max it - D(r,v))) -p m 6T max u, 



(Ps) 



(167) 



where p r , p v and p m are the costate variables associated to r, v and to, respec- 
tively. The switching function is 



$ := D U H = T max ((l -p m )b + Pv/m). 



(168) 
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Hence, the bang-bang optimal control is given by 

4=0 if $ f > 0, 
i* t =l if $ t < 0, 



(169) 



and the singular control can be obtained by formally solving $ = 0. The 
expression of u* ingular , however, is quite complicated and is not recalled here. 
The solution for this problem has the well-known typical structure 1-singular- 
0, as shown on Figures [3] and |4j 



ALTITUDE (NORMALIZED) 



VELOCITY (NORMALIZED) 





MASS (NORMALIZED) 




0.02 0.04 0.0 



Figure 3: Goddard Problem 

Shooting formulations. Once again fixing the control structure, the shooting 
unknowns are 

v = (pi,o s P2,o,P3,o,*i,*2,T) € M 6 . (170) 
Here is the classical shooting formulation with the entry conditions on t\ 

S 3 {u) := (x hT - 1.01, p 2>T ,p 3 ,T + l,$t 1 ,®t x ,H T ), (171) 

while the extended formulation is 

S 3 (y) := (xi, T - 1.01,pa,T,P3,T + 1, ® tl ,K>H T , [H] tl ,[H] t2 ). (172) 

10.2 Results 

All tests were run on a 12-core platform, with the parallelized (OPENMP) 
version of the SHOOT ([32)) package. The ODE solver is a fixed step 4th. 
order Runge Kutta method with 500 steps. The classical shooting is solved 
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SWITCHING FUNCTION 



6000r 
4000- 
2000 


-2000 
-4000 - 
-6000 - 
-8000- 
-10000- 
-12000' 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 
TIME 



Figure 4: Goddard Problem 



with a basic Newton method, and the extended shooting with a basic Gauss- 
Newton method. Both algorithms use a fixed step length of 1 and a maximum 
of 1000 iterations. In addition to the singular/bang structure, the value of the 
control on the bang arcs is also fixed according to the expected solution. 

The values for the initial costates are taken in [—10, 10], and the values for 
the entry/exit times in [0, T] for (Pi) and (P2). For (P3), the entry, exit and 
final times are taken in [0,0.2]. The number of grid points is set around to 10 
000 for the three problems. These grids for the starting points are quite large 
and rough, which explains the low success rate for (Pi) and (-P3). However, the 
solution was found for all three problems. 

For each problem, the results are summarized in 3 tables. The first table 
indicates the total CPU time for all shootings over the grid, the success rate of 
convergence to the solution, the norm of the shooting function at the solution, 
and the objective value. The second table recalls the solution found by both 
formulations: initial costate and junction times, as well as final time for (P3). 
The third table gives the singular values for the Jacobian matrix at the solution, 
as well as its condition number k :— a\/a n . 

We observe that for all three problems (Pi), (P2) and (P3), both formula- 
tions converge to the same solution, v* and the objective being identical to more 
than 6 digits. The success rate over the grid, total CPU time and norm of the 
shooting function at the solution are close for both formulations. Concerning 
the singular values and condition number of the Jacobian matrix, we note that 
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for (P2) the extended formulation has the smallest singular value going from 
10~ 8 to 1, thus improving the condition number by a factor 10 8 . This is caused 
by the combination of the two entry conditions into a single one that we used 
in the classical formulation for this problem: as the singular arc lasts until tf, 
there is only one additional unknown, the entry time. 

Overall, these results validate the extended shooting formulation, which per- 
form at least as well as the classical formulation and has a theoretical foundation. 

Remark 10.3. Several additional tests runs were made using the HYBRD 
( 14$ ) and NL2SNO solvers for the classical and extended shootings in- 

stead of the basic Newton and Gauss-Newton method. The results were similar, 
apart from a higher success rate for the HYBRD solver compared to NL2SNO. 

Remark 10.4. We also tested both formulations using the sign of the switching 
function to determine the control value over the bang arcs, instead of forcing the 
value. However, this causes a numerical instability at the exit of a singular arc, 
where the switching function is supposed to be but whose sign determines the 
control at the beginning of the following bang arc. This instability leads to much 
more erratic results for both shooting formulations, but with the same general 
tendencies. 

Problem 1: 

Shooting grid: [-10,10] x [0,T] 2 , 21 3 gridpoints, 9261 shootings. 



Shooting 


CPU Success Convergence Objective 


Classical 
Extended 


74 s 21.28 % 1.43E-16 -106.9059979 
86 s 22.52 % 6.51E-16 -106.9059979 



Table 1.1: (Pi) CPU times, success rate, convergence and objective 



Shooting 


Po h h 


Classical 
Extended 


-0.462254744307241 2.37041478456004 6.98877992494185 
-0.462254744307242 2.37041478456004 6.98877992494185 



Table 1.2: (Pi) solution v* found 



Shooting 


a\ er 2 03 


K 


Classical 
Extended 


3.61 0.43 5.63E-02 
27.2 1.71 3.53E-01 


64.12 
77.05 



Table 1.3: (Pi) singular values and condition number for the Jacobian 



Problem 2 

Shooting grid: [-10, 10] 2 x [0,T], 21 3 gridpoints, 9261 shootings. 
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Shooting 


CPU Success Convergence Objective 


Classical 
Extended 


468 s 94.14 % 1.17E-16 0.37699193037 
419 s 99.36 % 1.22E-13 0.37699193037 



Table 2.1: (P 2 ) CPU times, success rate, convergence and objective 
Shooting 



Pifl 



P2,0 



Classical 
Extended 



0.942173346483640 1.44191017584598 1.41376408762863 
0.942173346476773 1.44191017581021 1.41376408762893 



Table 2.2: (P 2 ) solution v* found 



Shooting 




K 


Classical 
Extended 


24.66 5.19 1.96E-08 
24.70 5.97 1.13 


1.26E+09 
21.86 



Table 2.3: (P 2 ) singular values and condition number for the Jacobian 



Problem 3 

Shooting grid: [-10, 10] 3 x [0,0.2] 3 , 4 3 x 5 3 gridpoints, 8000 shootings. 



Shooting 


CPU 


Success 


Convergence 


Objective 


Classical 


42 s 


0.82 % 


5.27E-13 


-0.634130666 


Extended 


52 s 


0.85 % 


1.29E-10 


-0.634130666 



Table 3.1: (P 3 ) CPU times, success rate, convergence and objective 



Shoot. 


Pr,0 


Pv.O 


Pm.O 


Class. 
Exten. 


-50.9280055899288 
-50.9280055901093 


-1.94115676279896 
-1.94115676280611 


-0.693270270795148 
-0.693270270787320 




h 




*/ 


Class. 
Exten. 


0.02350968417421373 
0.02350968417420884 


0.06684546924474312 
0.06684546924565564 


0.174129456729642 
0.174129456733106 


Table 3.2: (P 3 ) solution v* found 



Shooting 


cti 


ct 2 


ct 3 


0-4 


^5 


0-6 


K 


Classical 


6182 


9.44 


8.13 


2.46 


0.86 


1.09E-03 


5.67E+06 


Extended 


6189 


12.30 


8.23 


2.49 


0.86 


1.09E-03 


5.67E+06 



Table 3.3: (P3) singular values and condition number for the Jacobian 
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11 Conclusions 



Theorems 5.8 and 8.8 provide a theoretical support for an extension of the 
shooting algorithm for problems with all the control variables entering linearly 
and having singular arcs. The shooting functions here presented are not the 
ones usually implemented in numerical methods as we have already pointed 
out in previous section. They come from systems having more equations than 
unknowns in the general case, while before in practice only square systems have 
been used. Anyway, we are not able to prove the injectivity of the derivative of 
the shooting function when we remove some equations, i.e. we are not able to 
determine which equations are redundant, and we suspect that it can vary for 
different problems. 

The proposed algorithm was tested in three simple problems, where we com- 
pared its performance with the classical shooting method for square systems. 
The percentages of convergence are similar in both approaches, the singular 
values and condition number of the Jacobian matrix of the shooting function 
coincide in two problems, and are better for our formulation in one of the prob- 
lems. Summarizing, we can observe that the proposed method works as well as 
the one currently used in practice and has a theoretical foundation. 

In the bang-singular-bang case, as in the fishing and Goddard's problems, 
our formulation coincides with the algorithm proposed by Maurer [5]. 

Whenever the system can be reduced to a square one, given that the suffi- 
cient condition for the non-singularity of the Jacobian of the shooting function 
coincides with a sufficient condition for optimality, we could established the 
stability of the optimal local solution under small perturbations of the data. 
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